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ABSTRACT 


This report deals with progress made on the Grant NSG-3048 during the 
calendar year beginning March 1, 1978 and ending February 28, 1979. The 
NASA Technical Officer for this period was Dr. Kurt Seldner of Lewis Re- 
search Center. The directors of the research at the University of Notre 
Dame were Dr. Michael K. Sain, for the duration of the year and Dr. R. 
Jeffrey Leake, for the initial six month period. 

General goals of the research have been classified in two categories. 
The first category involves the use of modern multivariable frequency do- 
main methods for control of engine models in the neighborhood of a set- 
point. The second category involves the use of nonlinear modelling and 
optimization techniques for control of engine models over a more extensive 
part of the flight envelope. 

Progress in the first category has included the extension of CARDIAD 
(Complex Acceptability Region for DIAgonal Dominance) methods developed 
with the help of the grant to the case of engine models with four inputs 
and four outputs. A suitable bounding procedure for the dominance function 
has been determined. The bound is quadratic in the compensator elements, 
and can be helpful in the general case of multiple inputs and outputs even 
when the Hessian matrix is indefinite. In addition, improvements have 
been made on the 370 CARDIAD software; and a beginning in the process of 
developing interactive CARDIAD software on a recently acquired PDP-11/60 
has been made. These steps are expected to enhance grant impact in the 
area of noninteractive engine control design. 
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Progress in the second category has had its principal focus on auto- 
matic nonlinear model generation. The approach consists of using the gen- 
eral form 

x - A(x,u) (x - g(u)> 

in which steady state measurements of x, A, and B in a linearized model 
are used to determine values of A(x,u), g(u) , and 3g(u)/8u at setpoints. 
A Hermit e polynomial interpolation scheme has been employed to determine 
global extensions of A(x,u) and g(u) . Simulations of these models have 
produced satisfactory results where compared with the NASA DYNGEN digital 
engine deck. Studies have been begun to apply time-optimal control compu- 
tational techniques to the new models and to compare these results with 
those of prior studies under this support. 



I. GENERAL BACKGROUND 


The following sections describe researches under Grant NSG-3048 during 
the calendar year March 1, 1978 - February 28, 1979. In this section, we 
make a number of preliminary and general remarks and lay some of the ground 
work for those sections. 

Initiation of Grant NSG-3048 in March 1975 was timed with develop- 
ments in the engine industry, which was beginning to experience some lim- 
itations in the application of classical hydromechanical control technique 
as the primary base technology for modern engines with ever increasing 
sophistication. At the same time, milestone developments in digital 
hardware began to open realistic possibilities for onboard computation 
to an extent not heretofore possible. This confluence of events led di- 
rectly to the concept of increasing the role of electronics in engine 
control. In turn, the availability of digital electronics itself created 
a wide variety of opportunity for application of new control design phil- 
osophy and technique. Among the earliest of such studies is the F100 
Multivariable Control Synthesis Program [1] sponsored by the National 
Aeronautics and Space Administration, Lewis Research Center and the Air 
Force Aero-Propulsion Laboratory, Wright -Patterson Air Force Base. This 
program has recently completed the test phase. 

The advent of digital technology on the engine scene offers not 
only the opportunity to control more engine variables but also the pos- 
sibility of integrating engine and airframe control. Studies of this 
type have also begun. 
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Primary tools in the F100 Multivariable Control Synthesis Program 
were linear quadratic regulator (LQR) theory in the linear case. For the 
global control, nonlinear optimal methods were not directly applied. 

The purpose of Grant NSG-3048 has been to evaluate alternatives to 
LQR in the linear case and to examine nonlinear modelling and optimization 
approaches for global control. 

Context for the studies is set by the DYNGEN digital simulator [2], 
Based upon earlier computer codes GENENG [3] and GENENG II [4], DYNGEN 
has the combined capabilities of [3] and [4], for calculating steady- 
state performance, together with the further capability for calculating 
transient performance. DYNGEN uses a modified Euler method to solve the 
differential equations which model the dynamics of the engine. This mod- 
ified Euler method permits the user to specify large time steps, for ex- 
ample a tenth of a second; and this can result in considerable savings of 
execution time. On the other hand, convergence problems are sometimes 
encountered with DYNGEN when small time steps are used. 

The DYNGEN digital simulation is particularized to a given situation 
by a process of loading data for the various maps associated with a given 
engine. The maps for the Grant NSG-3048 have been provided by engineering 
personnel at Lewis Research Center. These maps correspond to a paper 
engine, which is not closely identified with any current engine. But 
the data do correspond in a broad, general sense to realistic two spool 
turbofan engines. The simulation provides for two essential controls, 
main burner fuel flow and jet exhaust area. Portions of the envelope 
which can be used for linear or nonlinear experimentation are limited by 
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the convergence capabilities of the available engine data on DYNGEN. 

As mentioned in the Final Report for NASA Grant NSG-3048, Supplement 
No. 1, a promising new technique for designing dynamical compensation be- 
gan to develop in the Fall of 1976. This methodology, built upon what are 
currently being called CARDIAD plots, was only being tentatively consid- 
ered in October, 1976 when the continuation proposal for NASA Grant NSG- 
3048, Supplement No. 2, was being written. Based upon favorable prelim- 
inary reaction by personnel from NASA Lewis Research Center, a decision 
was made to investigate further the use of CARDIAD plots as a design aid 
for turbofan engine control in the frequency domain. In essence, this 
study proved to be successful enough that it really dominated the re- 
maining time period of Supplement No. 1 and has continued through seccessive 
study periods. 

A great deal of the power of the CARDIAD plot arises from its sim- 
plicity. For each frequency, a circle is constructed on a planar plot. 

Data for the center and radius of this circle is obtained from the com- 
plex transfer function matrix of the plant. The circle may be solid or 
dashed. If solid, the inside of the circle defines the acceptable com- 
plex region for the value of a frequency dependent compensator element 
iti order to achieve dominance. If dashed, the outside of the circle de- 
fines the acceptable region. As the frequency follows a standard Nyquist 
pattern, these circles result in a CARDIAD plot. (Complex Acceptability 
JRegion for PI Ag onal Dominance) . This plot has been shown to speak con- 
structively to the issue of compensator choice to reduce interaction. 

Examples of the use of CARDIAD plots utilizing a two-input, five- 
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state, two-output engine model in which the inputs are fuel flow and noz- 
zle area, the states are compressor speed, fan speed, burner exit pres- 
sure, after-burner exit pressure, and high turbine inlet energy, and the 
outputs are thrust and high turbine inlet temperature, may be found in 
the Final Report for Grant NSG-3048, Supplement No. 2. 

Though linear multivariable frequency domain techniques are very prac 
tical as a part of an overall engine design, there is always the task of 
integrating them into an overall controller which is nonlinear and which 
operates over a substantial flight envelope. One approach to this need 
has been to design a type of model following control system in which a 
nonlinear model generates trajectories which can then be tracked by the 
actual engine. Linear multivariable controllers are often helpful for 
local adjustments to this tracking process. It follows that the model 
which generates the trajectories to be tracked is of major importance in 
any design. 

Some desirable features of such models are accurate steady-state be- 
havior and acceptable dynamical behavior. After consideration of several 
modelling ideas over the history of this study. Dr. R.J. Leake and J.G. 
Comiskey have suggested a form 

x = A(x,u) (x - g(u)) 

where g(u) can be developed to yield steady-state accuracy and A(x,u) 
can be determined for dynamical purposes. Moreover, the process of 
gathering data to identify these functions is well suited to digital eng- 


ine simulations 


II. SUMMARY OF RESEARCH RESULTS 

This section provides a brief technical summary of the research re- 
sults associated with the calendar year ending February 28, 1979. Further 
Information on the ideas discussed here may be found in Sections III and 
IV, and in their referenced appendices. 

1. The CARDIAD design method was applied to a turbofan engine model having 
four Inputs, four outputs, and six states. This was a major test for 

a graphical, interactive method, inasmuch as the engine model has only 
two major inputs, with the other inputs intended for a less major role. 
The application was successful. 

2. A suitable quadratic bound for the diagonal dominance function has been 
selected. This should make it possible to study the CARDIAD method in 
a general n- input, n-output case. 

3 . It has been found that the bound above can be governed in engine mod- 
els by a Hessian matrix which is indefinite. A method to design in 
that case has been proposed. 

4. Software for CARDIAD design on the 370 computer has been improved and 
modularized. The CARDIAD method is being prepared for use on a PDP-11/60 
computer with graphics interface. This is expected to reduce design 
time by orders of magnitude. Previous CARDIAD plots came by Calcomp 

on the 370 machine. 

5. The nonlinear modelling class 

x » A(x,u) (x - g(u) ) 

described in the previous Annual Technical Report has been applied to 
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the DYNGEN digital simulator. Responses of the models agree well with 
those of the digital deck. 

6. A Hermite interpolation scheme has been proposed to extend this class 
of models to a region of finite extent on the available DYNGEN envelope. 
Application of the method indicates acceptable results, but only with- 
in the region of the original data. 

7. A program of testing these models under time-optimal feedback control 
i and comparing their performance with that of the DYNGEN deck has been 

begun. 

8. Numerous technical improvements have been made in the mathematical pro- 
gramming for the time-optimal problem. These will be reported in the 
final version of an M.S. Thesis by J.G. Comiskey* 


I 


I 

I 


I 


I 

I 


I 
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III. LOCAL MULTIVARIABLE FREQUENCY DOMAIN METHODS 

In accordance with discussions carried out with engineering personnel 
at Lewis Research Center, a decision was made in the latter part of summer, 
1978, to begin a gradual phasing out of grant researches having to do with 
multivariable frequency domain methods. The motivation for this decision 
grew out of the feeling that grant activities had provided already a sub- 
stantial stimulus to work in this area. See, for example. Grant NSG-3048 
Annual Technical Reports for the last three years. Correspondingly, an in- 
crease of emphasis upon nonlinear modelling and control was stated as a 
grant goal. 

Because of the promise afforded by the CARDIAD method to achieve dia- 
gonal dominance in turbofan engine models, it was also decided to place 
the remaining frequency domain emphasis in the grant upon this method. 

Primary efforts have been those of R.M. Schafer, under the direction 
of Dr. M.K. Sain. The report on this work is divided into two six-month 
portions. 

1. First Six Months 

During the six-month period beginning on March 1, 1978, the CARDIAD 
method for design to achieve diagonal dominance has been advanced to the 
case in which plants have four inputs and four outputs. The following 
paragraphs review this progress and illustrate it with design examples. 

The figures for this discussion are many in numbers, and have been 
included as Appendix B. — 

Recall that in the CARDIAD plot approach, column dominance may be 
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achieved by using only precompensation. The general form of this precom- 
pensator, K(s) is given by 
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21 (s) 

1 
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k 32 (s) 
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k 34 

41 (S) 

k 42 (s) 

k 43 ( s) 
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An exception to this general form is the use of a K matrix to interchange 
columns of the system to facilitate achieving dominance. This procedure 
is ^explicated in the first design example. 


Let G(s) be the transfer function matrix of the uncompensated sys- 
tem. With the above form of K(s) , dominance in each column of the com- 
pensated system, Q(s) = G(s) K(s), is a function of the off-diagonal en- 
tries in the corresponding column of K(s) . Hence, the goal is to design 
the off-diagonal entries of a column of K(s) in such a way that Q(s) is 
dominant in that column. 


To choose the off-diagonal entries, the system is evaluated over a 
range of frequencies. At each frequency, the off-diagonal entries in each 
column of the compensator K(s)j g __. w are complex numbers; and a sufficient 
condition for dominance can be expressed as an inequality. If is de- 
fined as a vector of the real and imaginary parts of the off-diagonal 
entries in a column, then that column of the system will be dominant at 
the frequency of evaluation if 


f(£) > 0, 


where f(^c) is a function of the system evaluated at a frequency. ___ 

Two approaches are used to maximize f(^) and achieve dominance. The 
first approach is to draw the CARDIAD plots for each off-diagonal entry 


11 


in a column under the assumption that the other off-diagonal entries in 
that column of K(s) are zero. This is called the standard or Type 1 
analysis and is a partial gradient analysis of f(j^). Implicit in this 
type analysis is the fact that it is being attempted to achieve dom- 
inance using only one off-diagonal entry in each column. When dominance 
cannot be achieved using Type 1 analysis, a full gradient or Type 2 
analysis is performed on f(^). In Type 2 analysis, all off-diagonal 
entries are used to achieve dominance in the column. 

i 

In each type of analysis, the CARDIAD plots describe the choice of 
off-diagonal entries that will achieve dominance. Before proceeding witH 
the design examples, a review of the features of CARDIAD plots is in 
order. 

Consider first the Type 1 analysis. The extrema found in the CARDIAD 
plot analysis are plotted in the complex plane using one of three different 
symbols, two of which have circles associated with them. These are a 
*+' with a solid circle around it and an 'x' with a dashed circle around 
it. The former case occurs when a maximum is found and is such that 
f(x) > 0. The solid circle then encloses the region where the 

choice of any complex value inside the circle will achieve dominance at 
the frequency associated with that circle. Similarly, the 'x' and the 
dashed circle are drawn when the acceptable region lies outside the circle. 
Hence, in the Type 1 analysis, choosing an off-diagonal entry k^(s) 
which, at each frequency, is inside the circle if it is solid, or outside 
the circle if it is dashed, will achieve dominance. In the simplest case, 
there is a point on the real line which satisfies the above condition, 
and dominance can be achieved in the column with constant compensation. 
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It Is not always possible to achieve dominance using only one off- 
diagonal entry in a column. At frequencies at which this occurs, a 'A* 
is plotted. The 'A* indicates the best choice of a compensator at a 
given frequency, but the one entry alone will not be sufficient to make 
f(jc> > 0. 

The occurrence of this ’A’ symbol in Type 1 analysis plots will have 
one of 2 effects. If, for instance, some of the plots for a column of 
the system contain triangles, and others do not, dominance can be achieved 
by choosing an entry to fit one of the plots without triangles and let- 
ting the other entries remain zero. If, however, all the plots for a 
given column contain some triangles, no one entry will be sufficient to 
achieve dominance; and Type 2 analysis must be used. 

In the Type 2 analysis, the type of the center and circle drawn is 
decided in the same way as in Type 1 analysis, that is, under the assump- 
tion that the other off-diagonal entries are zero. The reason the grad- 
ient values are not used in this determination is to avoid misleading the 
designer. If the gradient values were used, the region indicated by a 
solid circle in one plot would only be valid of the centers of the other 
plots for that column were fit exactly. Hence, when using Type 2 analysis, 
the strategy is to fit the compensators to all the plots and then use 
Type 1 analysis to check for dominance. 

Notice carefully that, although the type (+,x,A) of the center and 
the type (solid, dashed) of the circle are determined by Type 1 analysis, 
the location of the center is determined by full gradient (Type 2) analysis. 


We now turn to the discussion of some design examples 


The model used for the design examples was taken from a paper by 

Peczkowski and Sain [6] . A state feedback model is used with Inputs 

WFMB, A., CIW, RCW; states N. , N„, P_, T. , T. c , . The states fed 
j 1 2* 7 4.5hi 4.51o 

back are N, , N_, P_, T. . and T, . This model was made dominant 
1* 2 7* 4.5hi 4.51o 

using CARDIAD plot analysis with two different approaches. 

In the first approach, the first step in the design procedure was 
to switch the first and third columns by choosing as a first compensator 

K « 

While this choice happens to make the first column dominant, it was chosen 
to facilitate achieving dominance in the third and fourth columns. Re- 
gardless of what column switch is used, dominance is easily achieved in 
the first two columns of most jet engine models since fuel flow and ex- 
haust area have a dominant effect on all states and outputs of a jet engine 
The improvement gained by the column switch, as will be seen later, is the 
simple shapes of the Type 2 analysis of the third column. 

The Type 1 analysis plots for the system with the column switch are 
given in Figures 1-24. Analysis of these plots proceeds as follows. 

The first column of the system, Figures 1 - 6, is dominant without 
further compensation. This can be seen from the fact that in each plot, 
the origin is included by all solid circles and excluded by all dashed 
circles. Thus, the plots indicate that an acceptable choice for any one 
off-diagonal entry (assuming the other off-diagonal entries are zero since 


0 0 10 
0 10 0 
10 0 0 
0 0 0 1 



this is Type 1 analysis) is zero. Hence, all further compensation of the 
first column will be the first column of the identity matrix. 

The second column. Figures 7-12, is not dominant; and the (1,2) 

' and (3,2) plots each have triangles at high frequencies, indicating that 
dominance cannot be achieved by using either one of these entries by 
themselves. However, the (4,2) entry. Figures 11 & 12, contains no tri- 
angles; and therefore dominance can be achieved by judicious selection of 
this entry. 

In the third column. Figures 13 - 18, all the plots contain triangles. 
This says that no one entry in the column will achieve dominance, and 
Type 2 analysis will be necessary. 

The fourth column plots. Figures 19 - 24, also all contain triangles; 
but the (2,4) entry has solid or dashed circles at frequencies of u _> 40, 
Hence, dominance can be achieved at these higher frequencies by designing 
for this entry. 

The Type 2 analysis plots of column 3 are given in Figures 25 - 30. 
Recall that when using Type 2 analysis, all three off-diagonal entries 
will be designed. The semi-circular shapes of the Type 2 analysis plots 
indicate that designing to fit the shapes will not be difficult. These 
easily designed shapes are the reason that the column switch was chosen 
as it was . 

From this analysis, the following entries in the compensator were de- 
signed. 

The (4,2) entry of K(s) was chosen to be 
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k 42 <s) 


,97s - 89.1 
.033s + 1 


Referring to Figures 11 and 12, this entry, as a function of frequency, 
starts on the real line to the left of the low frequency dashed circles, 
and moves clockwise through the complex plane such that at each. frequency , 
the value is outside the circle at that frequency if it is dashed, or in- 
side if it is solid. Note that there is no point on the real axis which 
is inside all solid circles and outside all dashed circles. Hence, there 
is no constant entry which will suffice. 

The three off-diagonal entries in the third column were chosen to 
be first order and designed to fit the shapes of the triangles of the 
Type 2 analysis plots (Figures 25 - 30) . The resulting entries are 

-.00019s - .0084 
k 13 (s) “ -.01137s + 1 

, t \ -.00008s - .00049 

*23 “ -.0113s + 1 

. f v .00003s + .000137 
43 ls; ~ -.0108s + 1 


Finally, the (2,4) entry was designed to achieve dominance for in >_ 40. 
This entry was designed to start at the w * 0 triangle, center on the 
hi E 40 solid circle, and remain outside the dashed circles which occur 
thereafter. The resulting entry is 

, , x -.00806s + .0212 

k 24 (s) * ""-Toms' +" i — • 


The resulting compensator is 


I^Cs) 


.97s - 89.1 
,033s + 1 


-.00019s - .0084 
-.01137s + 1 

-.00008s - .00049 
-,0113s + 1 


.00003s + .000137 
-.0108s + 1 


-.00806s + .0212 
-,0172s + 1 


0 


and the overall compensation thus far is K(s) * K K^(s). This compensator 
Is labeled KMP3 on the plots. 


The Type 1 analysis plots of the system compensated with K(s) equal 
to K K^(s) are given in Figures 31 - 54. The plots of the first three 
Columns show that these columns are now dominant at all frequencies. This 
is true since in each CARDIAD plot of the first three columns, the origin 
is included by all solid circles and excluded by all dashed circles. The 
fourth column is not dominant at all frequencies, but is dominant at fre- 
quencies of to _> 40 as expected. Thus, further compensation is necessary 
only in the fourth column. 

Since the Type 1 analysis plots of the fourth column. Figures 49 - 54, 
all contain triangles at lower frequencies, Type 2 analysis will be used 
to further compensate the column. The Type 2 analysis plots are given in 
Figures 55 - 60. 

Considering the (3,4) entry (Figures 59,60) first, we see that all of 
the centers are located at approximately 57.3 on the real axis. Hence, 
this value was chosen for the (3,4) entry. In a similar spirit, the (1,4) 
entry was chosen to be the value at w * 0, namely .9. The next compensator 


chosen was 
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1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

1 

0 


.9 

0 

57.3 

1 


The overall compensation is now K(s) ** K K^(s) K£. The only column 
that is affected by is the fourth column. By considering the relative 
magnitudes of the fourth column entries of K. (s) K 7 , it was decided that 
an ^ approximation could be made. The effect of this approximation is that 
the zeros in the (1,4) and (3.4) entries of K^(s) could be replaced by .9 
and 57.3 respectively with the result being approximately the same as the 


full product result. Thus, the fourth column of K^(s) became 


.9 

.00806s + .0212 
-.0172s +1 - 

57.3 

1 

The fourth column was replotted with K^(s) changed as above. The 
plots are labeled KMP4 and are shown in Figures 61 - 66. 


Analyzing these plots, we see that the (1,4) entry now contains no 
triangles. Though the column is not yet dominant, dominance can now be 
achieved by fitting a frequency dependent entry to the solid circles of 
the (1,4) entry plot. The entry which fits these circles is 



(s) 


.232s + .4104 


.1127s + 1 
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and the next compensator was 


1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

1 

0 


,232s + .4104 
,1127s + 1. 

0 

0 

1 


Again* this only changes the fourth column which becomes 


,232s + .4104 . „ 

,1127s + 1 + - 9 

00806s + .0212 
0172s + 1 

57.3 

1 


The overall compensation is K(s) = K K(s) where 


and 


K(s) * 



,00019s - .0084 .334s + 1.3104 


u 

-.01137s + 1 

,1127s + 1 


-,00008s “ .00049 

-.00806s + .0212 

i 

-,0113s + 1 

,0172s + 1 

0 

1 

57.3 

,97s - 89.1 

.00003s + .000137 


,033s + 1 

-,0108s + 1 

X 


Figures 67 - 72 are the Type 1 analysis plots of the fourth column 
with the above compensation and are labeled KMP5. Dominance has been 


ach? eved 
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Using the same model, a second approach was used to achieve dominance. 
The second approach proved considerably more direct than the first approach. 
The starting point of the approach is to choose G ^(0) as a first compensa- 
tor Instead of the column switch that was used in the first approach. In 
this case, dominance was easily achieved by using this approach- 

With regards to the labeling of the plots, the following should be 
noted. Analysis Type 11 is analysis Type 1 with G ^(0) as a first compen- 
sator. Type 12 is Type 2 with G ^(0) as a first compensator. 


Figures 73 - 96 are the CARDIAD plot Type 1 analyses compensated by 


G*" 1 (0) . 


G'^O) - 


Analysis of the columns is as follows. 

The first column (Figures 73 - 78) is not dominant, but there exists 
a point on the real line in the plot of the (4,1) entry. Figures 77 - 78, 
which is included by all the solid circles. Hence, dominance can be 
achieved in this column by choosing = .076. 

The second column (Figures 79 - 84) is not dominant, but the plot of 
the (1,2) entry contains no triangles. A simple first order compensator 
was chosen to fit the solid circles of the plot. The entry chosen was 

k 12 (s) " Tl78s S + 1 * 


The third column plots. Figures 85 -90, show that the column is not 
dominant and that Type 2 analysis will be necessary since all plots con- 
tain triangles. 

The fourth column plots. Figures 91 - 96, show that the fourth column 
is dominant without further compensation. 


Figures 97 - 102 are the Type 2 analysis plots of the third column 
of the system already compensated by G ^(0), The plots are all very reg- 
ular in shape, and first order compensators were fit to each plot. The 
three entries chosen were 

-6.627s 


(s) 


.0639s + 1 


. r x -3.806s 
k 23 ^ " .0639s + 1 


k 43 (s) * ,0640s + 1 * 


The overall compensator is 


K<s) 


0 

.076 


.Is 


178s + 1 
1 


0 

0 


-6.627s 
.0639s + 1 

-3.806s 
,0639s + 1 

1 

-1.182s 
.064s + 1 


0 

1 


The plots for the system with compensation K(s) are given in Fig- 
ures 103 - 126. As can be seen from the plots, the system is dominant. 

It should be mentioned that the second of these two design examples 
was completed in about thirty minutes terminal time, though final plots 
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had to be awaited in a longer time frame. Moreover, the four Input, 
four output case Is definitely nontrivial for current models of jet engines 

Accordingly, the results support additional efforts to complete the 
development and refinement of CAKDIAD methods. 

2. Second Six Months 

During the six -month period beginning on September 1, 1978, work on 
the CARDIAD method was primarily divided into two areas. The first area 
involved a theoretical study of the dominance equation used in the analy- 
sis and a general proof of the bound that is employed. The second area in- 
volved software generalization and initial development work on the inter- 
active graphics software package. 


Before proceeding with the discussion of the theoretical aspects, 
some preliminary development is in order. Recall that the basis of the 
CARDIAD method is the definition of diagonal dominance which states that 
the j^ column of an n x n matrix Z(s) will be dominant at a point 
s e C if and only if 


■jj w 


n 

> I 

i=l 


Zy(s) 


( 1 ) 


Let G(s) be the transfer function matrix of the uncompensated 
plant. Let K(s) be a precompensator of a form having l's on the main 
diagonal and general frequency-dependent entries off the main diagonal. 

If Q(s) = G(s) K(s), then the condition for the j*"* 1 column of the com- 
pensated plant to be dominant as s is 
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q ij (s) 


( 2 ) 


Recall Chat Che Initial step in Che development of Che dominance 
equation by the CARDIAD method was to square (2) . This gives the following 
condition for dominance. 


q J3 U> 


n 

I 

i«l 


1 2 


q ij (s) 


(3) 


The next step in the procedure was to replace the condition in (3) 
with a simpler formulation involving a bound that results in a sufficient 
condition for dominance. Using the bound, the column of Q(s) will 

be dominant at s if 


| 2 - (n-1) J 1 


qij(s) 


> 0 . 


(4) 


To prove the validity of this bound in the general case, it is nec- 
essary to show that 

M , M , 

M l a 2 _> [ l a.] 2 (5) 

i=l 1 i«l 1 


for all real a ^ and for M * 2,3,4... To do this, we shall first show 


that 


M 

l 

i-i 


a M-l M n a i'l «% 

1 + l Z Cot + a. 2 ) > [ l a ] 

1 i-1 j=i+l x 3 -*-■* 1 


M 

l 

i-1 


( 6 ) 


This is a consequence of the following calculations . 


(a ± - c»j) > 0 
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2 2 

■* + cij 2a^ ttj 

Thus* with the aid of this observation, it follows that 
M-l M „ M-l M 


l l [a ± + a . ] > l l 2a a . 

i-1 j-i+1 1 J i-1 j-i+1 x J 


(7) 


Adding the same term to each side of (7) preserves the inequality; and we 
get the relationship 


M 

l 

i-i 


M-l M a n *■* o 

! «i + I l + “j 2 ! > l “i 

1 i-1 j-i+1 1 J 1 


M 

l 

1=1 


M-l M 

+ I I 2ot i a i * 

i=l j-i+1 1 3 


( 8 ) 


However, the right-hand member summation expression is 


M , M-l M M - 

l a + l l 2a a =[ l a } 2 . 

i-1 1 i-1 j-i+1 3 i=l 


(9) 


With the completion of this preliminary step, it remains to show that 

M , M , M-l M . _ 

M I a i = I + l I [a. + a, 1 

i-1 1 i-1 1 i-1 j-i+1 1 3 


( 10 ) 


for all M > 2. The proof is by induction. Consider the M = 2 case. 


For M = 2, (10) becomes 

2 <> 2 <> 1 2 

l _ r .2. v v r 2 


2 I “ ± " I “i + I I ' f« ± + O 

j 1 M j 1 * j 1 * J 


i-1 


i-1 


i-1 j-2 


(ID 


which becomes 


. 2 - 2 2 2 2 2 
2 oiji + 2 a2 = ci^ + a 2 + a l + a 2 

2 2 
- 2 a 1 +2 a 2 ^ » 


and thus the relationship has been shown for M = 2. Next we will show 
that, if (10) is assumed true for M _> 2, it is also true for M + 1; 
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and the proof will be complete. Asstame 


M . M _ M-l M . , 

M l “i ■ l <*i + 5! I (a + a ) . 

i-1 1 i-1 1 i-1 j*=i+l 1 J 


Adding identical terms to each member yields 


H l o, 2 + M « 2 - I « t 2 + l l (c 2 + a 2 > 

i-1 i-1 1 i=.l j«i+l 1 J 


+ M a 


M+l * 


Again, adding another term results in 



Now, manipulating the equation, we get 



thus, if (10) is true for M ^ 2, it is true also for M+l. 
by induction, the bound is valid. 


Hence, 
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th 


column of Q(s) 


Thus, we have as a sufficient condition for the j 
to be dominant at a frequency that 

Iq^Oo )! 2 * (n “ 1) I ^ij 2 > °- 

±*i 

After considerable algebraic manipulation, a more useful expression 
of the dominance equation can be shown to be of the form 


f j + c > 0. 


In this equation. A, B, and c are respectively a 2n-2 by 2n-2 matrix, a 2n-2 vector, 

and a scalar; and each are functions of the plant transfer function matrix 

evaluated at a specific frequency. All entries are real, and the matrix A 

th 

is Hermit ian. The subscript j indicates that dominance in the j 
column is being considered. This subscript will henceforth be dropped, 

. with the understanding that the equation refers to dominance in only one 
column. The vector is a 2n-2 vector composed of the real and Imaginary 
parts of the off-diagonal entries of the compensator in the column being 
considered. Consider, for example, the first column. The associated col- 
umn of the compensator matrix K(s) is: 

KjCs) — [ 1 , ^21 ^ * k^i ( s) , . . . ] . 


If we let k^Cs) 


, then the vector is 


s*=jco 

^ a 21* 0 21* a 31* 6 3i»***J * 


By choosing ^ such that f(x) is positive, one is selecting the values 
of the off-diagonal entries of the compensator that achieve dominance at 
the frequency being considered. 
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Within the framework of the quadratic form, it is easy to demonstrate 
how the various types of analysis are performed. In the case of Type 1 or 
partial gradient analysis, all of the entries in are set to zero ex- 
, cept for one (a,0) pair; then the gradient of f(^c) is set to zero, and 
a solution for the non-zero entries (a,$) is found. In the case of Type 
2 analysis, the standard gradient is taken. This results in the expression 

= (A + A C ) x + b = 0. 

Slixce A Is symmetric, the solution for £ is: 

^ - 1/2 A^b. 

One significant change was made in Type 2 analysis, and that was in 
the manner in which the circle types and radii are chosen. In the past, 
the centers of the circles were the gradient value but the type and size 
• of the circle was decided in the same manner as in Type 1 analysis. Thus, 
the type and size of the circle were decided under the assumption that the 
other off-diagonal entries were zero. In Type 2 analysis this is not the 
case, and thus the designer was advised to ignore the circles when using 
Type 2 analysis and just try to fit a frequency dependent entry to the 
centers. 

In the present software, the circle types and radii are decided in 
the following manner. For a given center in one plot, a deviation of a 
specified percentage is made in a worst case direction from the gradient 
values (i.e. circle centers) in all the other plots for the column. 

Then, the center type and radius are decided for the remaining entry. 
Consider a deviation of 10%. If a given plot has a very large circle at 
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a frequency, this implies that if the other entries are designed to he with- 
in 10% of the centers at this frequency, there is a lot of freedom in the 
remaining entry; and the designer does not have to fit this entry very pre- 
cisely. If a very small circle results, the designer has to be very pre- 
cise. In the case where a triangle results, it means that no value will 
achieve dominance in this entry at the frequency being considered if the 
other entries are 10% or more away from their gradient values. With this 
new way of deciding the circle types and radii, the designer is afforded 
more information as to how closely each of the entries need to be fit. 

A third means of making f(^c) positive was studied during this period. 
In Type 2 analysis, if the Hessian, which in this case is the matrix 2A, is 
negative definite, then the solution found by the gradient analysis is a 
global maximum. If the Hessian is positive definite, the solution is a 
global minimum. However, in the case where A is indefinite, the solution 
found is a saddle point; and Type 2 analysis cannot be used to achieve dom- 
inance. 

If A is indefinite, a solution for ^ that will make f(x) posi- 
tive is a scalar times the eigenvector of A associated with one of the 
positive eigenvalues of A. Let k be a positive scalar. Then 

- k 2 + k e*i> + •=. 

Because ^ is an eigenvector associated with a positive eigenvalue of A, 
the quadratic term will be positive; and the function can be made arbi- 
trarily positive. Thus, in the case of an indefinite Hessian, the eigen- 
vectors serve as a guide to a solution that achieves dominance. This method 
of achieving dominance is still under study and is in the process of being 
Implemented. 
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A considerable amount of time during this six-month period was spent 
on software. The first area of work involved the generalization of the 
370 version of the CARDIAD software. The principal changes involved 
generalizing the subroutine which solves for the values in A, b, and c 
of the quadratic form. Previous to this time, a different subroutine was 
called depending on the number of inputs and outputs of the problem; and 
the expressions for the entries were quite complex. At present, one sub- 
routine is able to solve for these values for any size plant, with the 

i 

only limiting factor being the size of dimensioned arrays. Other changes 
on this software included modularizing certain activities, polishing some 
code that was quickly written in the original software, and adding soft- 
ware to perform the new Type 2 analysis of circle type and radii. 

A large amount of time was involved with the development and initial 
implementation of an interactive graphics version of the CARDIAD software 
on a PDP-11/60 recently acquired by the Department of Electrical Engineer- 
ing. Much of the work was preliminary in nature because there was lead 
time involved with getting acquainted with the features of a new machine. 
Current plans include a package whereby the designer can draw the CARDIAD 
plot on the graphics terminal and fit a design to the plot, by selecting 
the order and type of complex entry that he desires and picking design 
values from the displayed plot. The designer will then be able to draw 
the locus of the entry over the CARDIAD plot to determine if it is a good 
candidate for achieving dominance. Having chosen an entry, the designer 
can then redraw the CARDIAD plot with the compensator and determine if 
dominance has been achieved. The CARDIAD method is by nature an inter- 
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active tool, and the necessity of having to wait for Calcomp plots has 
been a great hindrance in the past. It is hoped the interactive software 
will greatly reduce the time necessary to achieve dominance. 


IV. GLOBAL NONLINEAR OPTIMAL METHODS 


As discussed in the preceding section, the future investigations on 
this grant will begin to decrease the emphasis on multivariable frequency 
domain methods and start a gradual increase of the effort placed on non- 
linear work. 

The present studies underway on the DYNGEN digital engine simulator 
are to be completed. 

1. First Six Months 

During this period, the primary efforts were those of J.G. Comiskey, 
under the direction of Dr. R.J. Leake. 

With respect to nonlinear modelling and optimization, the emphasis 
has been twofold: to develop good analytical nonlinear models of the jet 

engine and to use these models in conjunction with techniques of mathe- 
matical programming in order to develop advances to global control over 
significant reaches of the flight envelope. 

In general, there are several aspects to this part of the investi- 
gation. First, it is possible to conceive the basic differential equa- 
tions from fundamental principles. In this case, there are usually about 
sixteen nonlinear differential equations, as well as a large number of 
nonlinear static functions which serve as part of the coupling between 
the equations. These functions often have more than one argument. If 
the equations arise in this fashion, then there is a significant need 
to identify the parameters. This must normally be done from the DYNGEN 
digital simulation. Second, it is possible to assume a general form for 
the nonlinear differential equations in such a way that fundamental pr in- 
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cllples are not Ignored but that added emphasis is placed upon general 
mathematical form. If this general form is chosen according to a scheme 
designed to make maximum use of the type of data which is directly avail- 
able from the digital simulation, then a type of ’’automatic” nonlinear 
model generation becomes possible. Third, whether the first or second 
modelling procedure is employed, there is almost always a need to con- 
sider the problem of reducing the order of the models. Though order 
reduction can often be highly mathematical in nature, it is almost al- 

i 

ways the case that the reduced order model depends upon the scaling of 
the equations. As a result, the final reduced models often depend in a 
nontrivial way upon physical insight, as well as mathematical method. 

Work on this grant has focused especially upon the first and second 
aspects of the modelling problem, with a gradual specialization toward 
automatic model generation. 

Insofar as optimization is concerned, the stress has been placed 
upon time optimal control, and considerable effort has been invested in 
specialized programming methodology designed to take maximum advantage 
of the particular features of jet engine models. 

The main effort during this six-month period has been in the area 
of model development. Fundamentals of the approach have been described 
in the Final Report for Supplement No. 2. For completeness, however, 
they are briefly sketched here. 

These models are based upon the general form 

x « f(x,u) 


where 
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f(x,u) - A(x,u) (x-g(u)). 


the equation 

x - g(u) 

, governs the important steady state analysis. Moreover, if a program such 
as DYGABCD [5] is available, it is shown in the Final Report for Supplement 
No. 2 that information concerning 

A(x,u) , g(u) , 9g(u)/9u 

can then be numerically generated from DYNGEN. 

» 

It should be noted that interpolation on g(u) from its points and 
derivatives involves a sort of Hermite problem, which is discussed in the 
following section. 

A feature of the approach is its authenticity with regard to equili- 
brium values, state matrix A, and DC gain values. 


At this point an M.S. Thesis was outlined for the purposes of com- 
paring this class of models with earlier modelling ideas studied on the 
grant and of using it for time optimal control calculations and simula- 
tions . 


The Table of Contents as originally planned is attached to this re' 
port as Appendix C. 

2. Second Six Months 

During this period, the primary efforts were again those of J.G. 
Comiskey, but this time under the direction of Dr. M.K. Sain. 

Production of the M.S. Thesis planned according to Appendix C was 
begun. Appendix D contains Chapters I through V of this document. The 
contents closely follow the plan of Appendix C. 
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A discovery of this part of the work is that the accuracy of the num- 
erical models can degrade rapidly outside the range of their data points. 
This means that the designer must be sure to include steady state and dy- 
namic data from the entire area of state and control space that he or she 
wishes to explore. 

These results are not unusual in the theory of polynomic approxima- 
tion in general. However, in view of the conceptual challenges encountered 
in the interpolation steps, it may be that the indications are toward the 
study of different types of nonlinear modelling methods. 

The beginning of such a new type of modelling idea is scheduled for 
the continuation of this grant. In addition, the completion of the Comiskey 
work will appear as part of the semi-annual report associated with that con- 


tinuation 


V. SPECIAL INITIATIVES RELATED TO GRANT WORK 


One of the biggest related efforts to work on Grant NSG-3048 was the 
International Forum on Alternatives for Linear Multivariable Control , 
sponsored by the National Engineering Consortium in Chicago during October 
1977. Dr. M.K. Sain served as Program Chairman and devised with the help 
of an advisory committee a Theme Problem based upon turbofan engine control 
for use by participants in the Forum. 

Dr. T.F. Edgar, Program Chairman for the 1979 Joint Automatic Control 
Conference, expressed interest in a continuation of the Forum idea at the 
JACC. Dr. Sain contacted all the Forum attendees and asked them to indi- 
cate their interest in the project. As a result, three sessions entitled 
"Further Alternatives for Linear Multivariable Control" are planned for 
the 1979 JACC. 

Preliminary information on these sessions is included here. 
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1979 JACC 


FURTHER ALTERNATIVES FOR LINEAR MULTIVARIABLE CONTROL 

Chairman and Organizer: Michael K. Sain 

University of Notre Dame 


Session I 


1. Regulator Design for the F100 Turbofan Engine 

S. Eng ell and N. Munro 
University of Manchester 

2. Frequency Dependent Precompensation for Dominance in a Four Input/ 
Output Theme Problem Model 

R.M. Schafer and M.K. Sain 
University of Notre Dame 

3. On Hidden Stability Margins in Multivariable Control 

Z.V. Rekasius 
Northwestern University 

4. Stability and Homotopy II — - 

R. Saeks and J. Murray 
Texas Tech University 


Session II 

5. Design of a Turbofan Engine Controller Via Eigenvalue/Eigenvector 

Assignment: A New Sensitivity Formulation 

S.R. Liberty, R.A. Maynard, and R.R. Mielke 
Old Dominion University 

6. Quasi-Upper Triangular Decomposition Applied to the Linearized Control 
of a Turbofan Engine — Further Results 

W.E. Holley and W. Chung 
Oregon State University 
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7. Design of Flight Control Systems Via Robust Decoupled Servomechanism 
Theory 

S.-H. Wang, University of Maryland and 
E.J. Davison, University of Toronto 

8. Computer Aided Design of Control Systems Via Optimization 

D.Q. Mayne 
Imperial College 


Session III 

9. Inverse Systems in Multivariable Control 

J.L. Peczkowski, Bendix Corporation, 

M.K. Sain, University of Notre Dame, and 
R.J. Leake, Fresno State University 

10. Failure Accommodation in Gas Turbine Engines 

R.K. Sahgal and R.J. Miller 
Pratt & Whitney Aircraft 

11. Model Algorithmic Control: Extensions and Further" Applications 

R.K. Mehra, A. Rault, and R. Rouhani 
Scientific Systems, Inc. 

12. An Application of Model Following Control 

J.D. Aplevich 
University of Waterloo 


VI. REFERENCES 



1. R.L. De Hoff, W. Earl Hall, Jr., R.J. Adams, and N.K. Gupta, "F100 
Multivariable Control Synthesis Program*', Vols. I and II, Technical 
Report AFAPL-TR-77-35, June 1977. 

2. J.F. Sellers and C.J. Daniele, "DYNGEN-A Program for Calculating Steady- 
State and Transient Performance of Turbojet and Turbofan Engines", NASA 
Technical Note TN D-7901, April 1975. 

3. R.W. Koenig and L.H. Fishbach, "GENENG-A Program for Calculating De- 
sign and Off-Design Performance for Turbojet and Turbofan Engines", 

NASA Technical Note TN D-6552, February 1972. 

4. L.H. Fishbach and R.W. Koenig, "GENENG II-A Program for Calculating 
Design and Off-Design Performance of Two and Three Spool Turbofans 
with as Many as Three Nozzles", NASA Technical Note TN D-6553, Feb- 
ruary 1972. 

5. "DYGABCD - A Program for Calculating A, B, C, and D Matrices for Linear 
Descriptions of Simulated Turbojet and Turbofan Engines", NASA Lewis 
Research Center, Preliminary Report. 

6. Joseph L. Peczkowski and Michael K. Sain, "Linear Multivariable Syn- 
thesis with Transfer Functions," in Alternatives for Linear Multivar- 
iable Control , M.K. Sain, J.L. Peczkowski, and J.L. Melsa, eds. Chicago: 
National Engineering Consortium, 1978, pp. 71-87. 



APPENDIX A 


GRANT BIBLIOGRAPHY, INCEPTION TO PRESENT 


(1) R.J. Leake, J.L. Melsa, and M.K. Sain, "Preliminary Proposal on Mod- 
ern Approaches to Jet Engine Control", November 1, 1974. 

(2) R.J. Leake, M.K. Sain, and J.L. Melsa, "Proposal to NASA for Support 
of a Work Entitled ’Alternatives for Jet Engine Control’", November 
27, 1974. 

(3) R.J. Leake and M.K. Sain, "Semi-Annual Status Report on NASA Grant 
NSG-3048, 'Alternatives for Jet Engine Control”', March 1, 1975 - 
August 31, 1975. 

(A) R.J. Leake, M.K. Sain, and J.L. Melsa, "Proposal for Continuation 
of NASA Grant NSG-3048, 'Alternatives for Jet Engine Control'", Oc- 
tober 1, 1975. 

(5) J.C. Shearer, "An IBM 370/158 Installation and User's Guide for the 
DYNGEN Jet Engine Simulator", M.S. Thesis, November 1975. 

(6) T.C. Brennan and R.J. Leake, "Simplified Simulation Models for Con- 
trol Studies of Turbojet Engines", Technical Report EE 757, Univer- 
sity of Notre Dame, November 1975. 

(7) R.J. Leake, M.K. Sain, and J.L. Melsa, "Final Technical Report on 
NASA Grant NSG-3048, 'Alternatives for Jet Engine Control'", March 1, 
1975 - February 29, 1976. 

(8) R.R. Gejji and M.K. Sain, "Polynomial Techniques Applied to Multi- 
variable Control of Jet Engines", Technical Report EE 761, March 1976. 

(9) V. Seshadri and M.K. Sain, "Multi variable System Compensation In- 
cluding Characteristic Methods for Jet Engine Control", Technical 
Report EE 763, May 1976. 

(10) A.J. Maloney III, "Graphics Analysis of Dominance in Jet Engine Con- 
trol Models", Technical Report EE 765, June 1976. 

(11) M.K. Sain, R.J. Leake, R. Basso, R. Gejji, A. Maloney, and V. Seshadri, 
"Alternative Methods for the Design of Jet Engine Control Systems", 
Proceedings Fifteenth Joint Automatic Control Conference , pp. 133- 
142, July 1976. 

(12) J.C. Shearer, W.E. Longenbaker, Jr., and R.J. Leake, "An IBM 370/158 
Installation and User's Guide to the DYNGEN Jet Engine Simulator", 
Technical Report EE 766, July 1976. 


38 


39 


(13) R.R. Gejji and M.K. Sain, "A Jet Engine Control Problem for Eval- 
uating Minimal Design Software", Proceedings Midwest Symposium on 
Circuits and Systems , pp. 238-243, August 1976. 

(14) V. Seshadri and M.K. Sain, "Interaction Studies on a Jet Engine 
Model by Characteristic Methodologies", Proceedings Midwest Sympo- 
sium on Circuits and Systems , pp. 232-237, August 1976. 

(15) R.J. Leake and M.K. Sain, "Semi-Annual Status Report on NASA Grant 
NSG-3048, ’Alternatives for Jet Engine Control', Supplement No. 1", 
March 1, 1976-August 31, 1976. 

(16) R. Basso and R.J. Leake, "Computational Alternatives to Obtain Time- 
Optimal Jet Engine Control", Technical Report EE 767, September 1976. 

(17) R. Basso and R.J. Leake, "Computational Methods to Obtain Time-Op- 
timal Jet Engine Control", Proceedings Fourteenth Allerton Conference 
on Circuit and System Theory , pp. 652-661, September 1976. 

(18) M.K. Sain and V. Seshadri, "A Result on Pole Assignment by Exterior 
Algebra", Proceedings Fourteenth Allerton Conference on Circuit and 
System Theory , pp. 399-407, September 1976. 

(19) R.R. Gejji, "A Computer Program to Find the Kernel of a Polynomial 
Operator", Proceedings Fourteenth Allerton Conference on Circuit 
and System Theory , pp. 1091-1100, September 1976. 

(20) R.J. Leake and M.K. Sain, "Proposal for Continuation of NASA Grant 
NSG-3048, 'Alternatives for Jet Engine Control'", October 25, 1976. 

(21) R.M. Schafer, W.E. Longenbaker, and M.K. Sain, "System Dominance: 

A Preliminary Report on an Alternate Approach", Technical Report, 
University of Notre Dame, November 25, 1976. 

(22) R.J. Leake and M.K. Sain, "Final Technical Report, NASA Grant NSG- 
3048, 'Alternatives for Jet Engine Control', Supplement No. 1", 

March 1, 1976 - February 28, 1977. 

(23) R.J. Leake, J.G. Allen, and R.S. Schlunt, "Optimal Regulators and 
Follow-Up Systems Determined from Input-Output Signal Monitoring", 
Proceedings IEEE International Symposium on Circuits and Systems , 
April 1977. 

(24) R.M. Schafer, "A Graphical Approach to System Dominance", Technical 
Report EE 772, University of Notre Dame, April 1, 1977. 

(25) W.E. Longenbaker and R.J. Leake, "Hierarchy of Simulation Models 
for a Turbofan Gas Engine", Proceedings Eighth Annual Pittsburgh 
Conference on Modeling and Simulation . April 1977. 

(26) P.W. Hoppner, "The Direct Approach to Compensation of Multivariable 


40 


Jet Engine Models", Technical Report EE 774, University of Notre Dame, 
May 1977. 

(27) R.R. Gejji and R.J. Leake, "Time-Optimal Control of a Single Spool 
Turbojet Engine Using a Linear Affine Model", Technical Report EE 
7711, University of Notre Dame, June 1977. 

(28) R.M. Schafer, R.R. Gejji, P.W. Hoppner, W.E. Longenbaker, and M.K. 
Sain, ’’Frequency Domain Compensation of a DYNGEN Turbofan Engine 
Model", Proceedings Sixteenth Joint Automatic Control Conference , 
pp. 1013-1018, June 1977. 

(29) R.R. Gejji and M.K. Sain, "Application of Polynomial Techniques to 
Multivariable Control of Jet Engines", Proceedings Fourth IFAC Sym- 
posium on Multivariable Technological Systems , pp. 421-429, July 
1977. 

(30) R.R. Gejji, "Use of DYGABCD Program at Off-Design Points", Technical 
Report EE 7703, University of Notre Dame, July 1977. 

(31) E.A. Sheridan and R.J. Leake, "Non-Interactive State Request Jet 
Engine Control with Non-Singular B Matrix", Proceedings Twentieth 
Midwest Symposium on Circuits and Systems , pp. 539-543, August 

1977. 

(32) R. Gejji, R.M. Schafer, M.K. Sain, and P, Hoppner, "A Comparison of 
Frequency Domain Techniques for Jet Engine Control System Design", 
Proceedings Twentieth Midwest Symposium on Circuits and Systems , 
pp. 680-685, August 1977. 

(33) W.E. Longenbaker and R.J. Leake, "Time Optimal Control of a Two-Spool 
Turbofan Jet Engine", Technical Report EE 7714, University of Notre 
Dame, September 1977. 

(34) R.J. Leake and J.G. Comiskey, "A Direct Method for Obtaining Non- 
linear Analytical Models of a Jet Engine", Proceedings International 
Forum on Alternatives for Linear Multivariable Control , National 
Electronics Conference, Chicago, pp. 203-212, October 1977. 

(35) M.K. Sain and R.J. Leake, "Proposal for Continuation of NASA Grant 
NSG-3048, ’Alternatives for Jet Engine Control*", October 25, 1977. 

(36) J.A. Ortega and R.J. Leake, "Discrete Maximum Principle with State 
Constrained Control", SIAM Journal on Control and Optimization , Vol. 
15, No. 6, pp. 109-115, November 1977. 

(37) Michael K. Sain and V. Seshadri, "Pole Assignment and a Theorem from 
Exterior Algebra", Proceedings IEEE Conference on Decision and Con- 
trol , pp. 291-295, December 1977. 

(38) R. Michael Schafer and Michael K. Sain, "Some Features of CARDIAD 
Plots for System Dominance", Proceedings IEEE Conference on Decision 


41 


and Control , pp. 801-806, December 1977. 

(39) M.K. Sain, "The Theme Problem", in Alternatives for Linear Multivar- 
iable Control , M.K. Sain, J.L. Peczkowski and J.L. Melsa, Editors. 
Chicago: National Engineering Consortium, 1978, pp. 20-30. 

(40) R.M. Schafer and M.K. Sain, "Input Compensation for Dominance of 
Turbofan Models", in Alternatives for Linear Multivariable Control , 
M.K. Sain, J.L. Peczkowski, and J.L. Melsa, Editors. Chicago: Na- 
tional Engineering Consortium, 1978, pp. 156-169. 

(41) J.L. Peczkowski and M.K. Sain, "Linear Multivariable Synthesis with 
Transfer Functions", in Alternatives for Linear Multivariable Con- 
trol , M.K. Sain, J.L. Peczkowski, and J.L. Melsa, Editors. Chicago: 
National Engineering Consortium, 1978, pp. 71-87. 

(42) R.J. Leake and M.K. Sain, "Semi-Annual Status Report, NASA Grant 
NSG-3048, 'Alternatives for Jet Engine Control', Supplement No. 2", 
March 1, 1977 - August 31, 1977. 

(43) R.J. Leake and M.K. Sain, "Final Technical Report, NASA Grant NSG- 
3048, 'Alternatives for Jet Engine Control', Supplement No. 2", 

March 1, 1977 - February 28, 1978. 

(44) M.K. Sain and R.J. Leake, "Semi-Annual Status Report, NASA Grant NSG- 

3048, 'Alternatives for Jet Engine Control'", March 1, 1978 - August 
31, 1978. — 

(45) R.J. Leake, J.L. Peczkowski, and M.K. Sain, "Step Trackable Linear 
Multivariable Plants", Technical Report EE 789, September 1978. 

(46) R.M. Schafer and M.K. Sain, "CARDIAD Design: Progress in the Four 

Input/Output Case", Proceedings Sixteenth Allerton Conference on 
Communication, Control, and Computing , p. 567, October 1978. 

(47) M.K. Sain, "Proposal for Continuation of NASA Grant NSG-3048, 'Al- 
ternatives for Jet Engine Control'", November 1, 1978. 

(48) M.K. Sain, "On Exosubsets and Internal Models", Proceedings IEEE 
Conference on Decision and Control, pp. 1069-1073, January 1979. 


(49) M.K. Sain, "Annual Technical Report, NASA Grant NSG-3048, 'Alter- 
natives for Jet Engine Control'", March 1, 1978 - February 28, 1979. 


Figure 1 


’t. ' 

' A'iV-- 

------ 


y J 

' ^ 




„ _ _ _ _ _ N_ 

-T — — tt v. 
s '«.'«,' V v 


-~'-~-~A 


R?> 


,f AOy- ^ \\\\ 

*\\ “ v “ cn> » X\ 

toi ' \ v v. -~ _. 

% V 1-N " „ 

\ s . ' _ 
|M$S ' 

W*\ >. ' - - 

V^nV 


■ V *. v. s ^ s N 

•is V \S V S S \ V 
V >>.\ ' s v O \ s 

v>v \ N V V V N 

\\ V W V S 

\\v \ vK\\ \ 

s\‘ #\\\\< 

H M / n <• \ i 

\ l I v V Wv 

t » l / ' w 

i / 1 • 1 y t 

\ /I' 1 • /' t 

•i y J « « 

\ * * a > < 1 

r l ~- ;/ - y ' " ; * * -+ < 1 

A'r-~ J' • I J' 

'/ y- "/ j~r / 

\-y<A > -''ii ‘ . 

]-'- /- III 1 

' - ' l/l 1 

*-< " ~ -* i. _ l/l 1 

" ~ y-‘ 

<'*■*» 'll/ j 

'' "'•-,'''// 

- S / / / ' 

' > -~r ! ! 

' * s " <_ ' 

' ! 

^ 


•0.03 -0.05 


0.07 o'. 11 oTTF 

REAL PART 


a»/a/v4t 


CARDIAC PLOT; ^ . 1 ENTRY: COLUMN SWITCH, 6/14/78. 


0F P0 °* «S£ INPUTS « 4, SYSTEM OROER * 6. ANALYSIS TYPE = 1 



0.09 -0.05 -0.01 0.03 0.07 0.11 0.1S 0.19 0.23 

REAL PART 

LOCUS OF CENTERS PLOT: 2.1 ENTRY: COLUMN SWITCH. 6/14/78 
INPUTS = 4. SYSTEM ORDER = 5. ANALYSIS TYPE = 1 






110.23 


-95.33 -91. M3 


-S7.03 -52.63 -33.23 -23.33 -9.H3 4.97 

REAL PRRT 


LOCUS OF CENTERS PLOT: 3,1 ENTRY: COLUMN SWITCH. 6/14/79 


INPUTS = 4. SYSTEM ORDER = 6. ANALYSIS TYPE = 1 























52 








-56. 3C -42.50 -23.40 -14.20 C.00 14.20 23.40 42.50 

REAL PART 


LOCUS OF CENTERS PLOT: 4.2 ENTRY : COLUMN 5WITCH. 6/14/78 


INPUTS = 4, SYSTEM ORDER = 6, ANALYSIS TYPE * 1 




54 


Figure 13 
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Figure 29 
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CARD I AO PLOT: 4.3 ENTRY: COLUMN SWITCH. 6/14/78. 
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Figure 31 
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Figure 35 
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Figure 45 
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Figure 50 
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Figure 56 
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Figure 59 
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Figure 61 








103 


Figure 62 



LOCUS OF CENTERS PLOT: 1,4 ENTRY; KMP 4, S/22/78. 
INPUTS = 4. SYSTEM ORDER = 6. ANALYSIS TYPE * I 


«■ 








REAL PART 



CARDIAO PLOT: 2,4 ENTRY: KMP 4, 6/22/79 
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Figure 64 
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LOCUS OF CENTERS PLOT; 2.4 ENTRY: KMP 5. 5/24/79. 
INPUTS = 4. SYSTEM OROER = 6. ANALYSIS TYPE = 1 
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Figure 72 





LOCUS OF CENTERS PLOT: 3,4 ENTRY : KMP 5, 6/24/78. 
INPUTS = 4. SYSTEM ORDER = 6. ANALYSIS TYPE = 1 
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CARO I AD PLOT: 3.3 ENTRT : 1 /G (0) . 6/24/78. 
INPUTS = 4, SYSTEM 0R0ER = 6. ANALYSIS TTPE = 11 
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INPUTS 4, SYSTEM ORDER = 6. RNflLYSIS TYPE = 11 
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Figure 107 
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Figure 110 
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CflRDIflD PLOT; 1,3 ENTRY: KMP l, 1/GCO), S/25/78. 
INPUTS -4, SYSTEM ORDER = 6. flNRLYSIS TYPE = 11 



157 




Figure 116 
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Figure 118 
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Figure 119 
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CHAPTER I 


INTRODUCTION 

The scope of this work will be to make preliminary efforts to generate 
non-linear numerical models of a two-spooled turbofan jet engine, and sub- 
ject these models to a known method of generating global, non-linear, time 
optimal control laws. The models will be derived numerically, directly from 
empirical data, as a first step in developing an automatic modelling proce- 
dure. 

A hierarchy of models, including analytical and numerical models, will 
be established. The numerical models will be described in detail, and their 
step responses compared to those of the hypothetical jet engine from which 
they were derived. A method of generating time optimal control laws will be 
explained, programmed, and applied to the numerical models. Finally, these 
control laws will be tested, both on the models from which they were gener- 
ated, and on the hypothetical jet engine. 

This is the third in a series of similar works, whose ultimate goal 
is the development of an automated modelling method. Even though DYNGEN, 
an elaborate mathematical model, already exists, new models were developed 
for two reasons. First, DYNGEN uses too much cpu time to be called re- 
peatedly by an iterative method such as dynamic programming; a smaller, 
faster model is required. Second, DYNGEN assumes the role of a physical 
plant in this work, since access to a real jet engine is impossible. 

In his paper. Basso [9] uses two methods of generating optimal control 
sequences. The first is the dynamic programming successive approximations 
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technique. This actually generates a control law, from which a control se- 
quence can be derived. The second is a modified FI etcher- Reeves conjugate 
gradient method. This method generates a control sequence that drives the 
system to the target in minimum time. The modification consists of the in- 
troduction of constraints into the original method. 

His findings were that both methods yielded similar results, and that 
tiie number of computations necessary to solve a problem increases geometric- 
ally with system order for dynamic programming, but only arithmetically for 
the conjugate gradient method. 

Longenbaker [1] applies the dynamic programming method to several models 
of the F-100 engine. His models include several linear systems, and one 
non-linear, analytical system of differential equations derived from phy- 
sical and mathematical relationships among the state and control variables. 
Longenbaker concludes that the agreement between this analytical model and 
the DYNGEN simulator is not strong enough to justify great faith in the 
control law generated. 

In this paper, the same dynamic programming method is applied with a 
modification introduced to reduce cpu time, to several numerical, non- 
linear models of the F-100 jet engine. The conclusions are that better 
numerical agreement with DYNGEN is achieved by this numerical models than 
l»y Longenbaker' s analytical model, with a much smaller expenditure of 
man hours. However, complex interpolation techniques cause these models 
to use extravagant amounts of cpu time. Either larger data bases and less 
interpolation, or a more economical technique like Basso's conjugate gra- 
dient method should be explored in future work. 


CHAPTER II 


TWO SPOOL TURBOFAN JET ENGINE MODELS 

2.1 Introduction 

A form was chosen for the system model, which isolates static and dy- 
namic portions of the system behavior, so that each of these can be modelled 
independently. Several methods of modelling these two portions were tried, 
resulting in a hierarchy of models. Each model was subjected to the same 
analysis for purposes of comparison. 

In this chapter, the system model form is derived, and the modelling 
methods outlined. These methods will be treated in detail in later chapters. 


2.2 Basic Approximation Approach [2] 

Now consider a method for obtaining nonlinear models. Let 


x = f(x,u) 


( 2 . 2 - 1 ) 


with x an n vector and u an m vector denoting a dynamical system such 
as a jet engine, in which the state variables and parameters u remain pos- 
itive throughout the system operation and there is a function g(u) such 
that for each equilibrium point 


f(x,u) = 0 < ► x = g(u) 


( 2 . 2 - 2 ) 


The steady state system analysis involves the study of the function g(u). 


We propose to approximate the system (2.2-1) by 


x=A(x,u)[x - g(u)] 


(2.2-3) 


where A(x) is a square matrix which varies as a function of x. Notice 
that if Xp is an equilibrium point of (2.2-1), 


*D " 8<V' 


then a lin- 
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earizatlon about this equilibrium point results in the linear system 

6x ® A^iSx + (2.2-?4) 

and a linearization of the approximating system (2) at x^ * g(u^) results in 

6x ® ACXjjJSx + [-ACXjP (u d )]6u (2.2-5) 

Hence, the linearization of (2) will match the linearization of (1) if and 
only if 

A< V ■ A d • -A D ^ %> ' B D (2 - 2 - 6) 

I 

Also, if A^ is invertible, as is often the case for jet engine models, 
equation yields 

|S (U D > = -P^\ (2.2-7) 

These static and dynamic data are available from known algorithms [7], 
leaving only the choice of interpolation methods for generating non-linear 
models . 

2.3 Hierarchy of Models 

This work has resulted in the formation of a hierarchy of models, each 
a step in the development of an automated modelling method. They are clas- 
sified as follows: 

Model 0: The actual F-100 type engine (hypothetical) 

Model 1: The DYNGEN [6] simulation program, coded with data presumed to 

have been taken from experimental measurements on Model 0. This 
model solves 16 nonlinear differential equations and uses data 
maps and thermodynamic tables which cannot be expressed analy- 
tically. 
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Model 2: The Longenbaker [1] model, a 5th order, nonlinear, analytical mod- 

el. It includes the 5 state differential equations which govern 
the dynamical behavior of the system, along with 20 algebraic 
equations which express the relationship between various engine 
variables. This model is discussed in detail in [1], 

Model 3A: The linear affine power law model, which is a fit of steady state 
data to a selected form with linear, nonlinear and constant terms. 

Model 3B: The straight linear affine model, generated in the same manner 
as 3A» without the non-linear terms, to serve as a comparison. 

Model 4: The Quasi-Hermite interpolation model. Also a fit to steady 

state data, this model employs value and derivative matching 
over a two dimensional subset of the state space. 

Models 3 and 4 will be outlined briefly here, ancTdetailed in later 
chapters . 

2.4 Linear Affine Power Law Model [2] 

This model approximates the system by interpolating values. of A(x,u) 

from values of the matrix at two data points, and by generating values, for 
g(u) by a fit of the form: 

C 3i C Ai 

g ± (u) = c^Uj + c 2l u 2 + c 5± u 1 u 2 + C 61 » 1 - 1, ... ,5 (2.4-1) 

to the same two data points. 

2.5 Quasi-Hermite Interpolation Model “ 

This model approximates the system by interpolating values of A(x,u) 
from values of the matrix at five data points, and by interpolating values 
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CHAPTER III 


LINEAR AFFINE POWER LAW MODEL 3 


3.1 Formation of A(x,u) [2] 

Values of A(x,u) are interpolated from the values of A^ * A(Xp,Up) 
and Ay *= A(Xy,Uy) in the following manner: 

A(x,u) = Ay d iag + A^ diag(^~ — (3.1-1) 


*Dj “*WJ 


D j *** 

th 


where diag(*) is a diagonal matrix which causes the j column of A(x) to 
be interpolated linearly between the jth columns of Ay and A^ with x^ 
as the interpolation variable. 


3.2 Approximation of g (u) [2] 

The parameter vector u is presumed to be made up of physical control 
variables, and parameters such as fuel flow and nozzle- area. The equilibrium 
function is to be approximated in a manner such that both the equilibrium 
values and the linearizations of the approximating system (3) match those 
of system (1) at both Xjy and x^. This requires then that 

* (u D> " *D 8< V = *W (3.2-1) 


and also 


£&(„)= _a -1 r i£ - * _1 ’ 

<5u 


( V ■ -*d b d • iSV'-VSw 


(3.2-2) 


The method we propose here is to approximate each scalar component g,.(u) 
of g(u) by a linear affine power law form 

g(u) - +. . .+ c n u m + o 2Bfl u 1 * fl u 2 m+2 . . . + c 2nif2 (3.2-3) 

for which the jth partial derivative is 
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SR i 

Su, 


C j + C 2m+l C m+j 


m+1 nrhn 

u ...u 
1 m 

u. 


(3.2-4) 


Now, If the variables are normalized and scaled such that 


Up ” (1,1, ...,1) = 1 u y = (a,a,...,a) = £ 
then, the conditions of (11) and (12) can be put in the form 

k j = “ c j + c 2m+l c m+j 

* k . . *= (a) = c + c, -c a 

trt+j 6Uj — j 2nrfl m+j 


k 2nrt-l S i^ " ^ C j + C 2nrfl + C 


2m+2 


J.< 


k 2m4-2 " 8 i (a) = a ^ c j + c 2nrH a ’ + c 


2m+2 


(3.2-5) 


(3.2-6) 


and summing the first two of these over j yields 


+ c 2m+l 


2 k mtJ = + c 2nri-l a 


£c.,-l 

h 




k 2mfl j + c 2m+l + C 2nri-2 


'h* + 


I< 


C rt » 1 ^ 


m+j 


+ c. 


2m+2 "+ w j 2m+l "2m+2 


(3.2-7) 


which is of the form 


S 1 = r l + r 3 r 2 


s 2 = r l + r 3 r 2 a 


V 1 


S 3 * r l + r 3 + r 4 


s. * ar, + r 0 a + r. 

4 13 4 


(3.2-8) 
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which. Incidentally Is the m=l condition also. This set of transcendental 
equations is solved numerically for r^, r 2 » r^, r^ and (3.2-6) is then 
used to solve for each c^. In the event that (3.2-8) has no solution, a 
best fit is made on the second equation by varying while the other con- 
ditions are satisfied exactly. 


3.3 Computational Algorithm [2] 

In this section, we present an algorithm which serves to automate the 
process of finding a nonlinear model for a system 


x * f(x,u) 


(3.3-1) 


to be approximated from Xp.Up.x^jU^ApjBj^A^B^ by a normalized system. 

The algorithm will automatically perform the normalization and, hence, ac- 
tually approximate the system 

£ = f(£,fi) ~~ (3.3-2) 


where £ 


i = x^/x^ » * u j/ u p • The approximating system is of the form 


& » A(£)[x-g(u)] 


where 


*i -*w 


iw ■ ^ Mag V7 + ^ dlag 


i i 


i 

i "i 


(3.3-3) 


(3.3-4) 


and 


ft. 


I C j U J + C^ n Tr ur j + c 


*m+j 


'2m+l . “j 


i 

2m+2 


where 


U * - oft + 


(3.3-5) 


(3.3-6) 


Algorithm I . 

1. Input: x D ,u D ,A D ,B D ,m,n,a,E,x w ,5i w ,^,B w 

where m = number of controls 


n ■ number of states 


179 


II 


I: 


2. Calculate: 

A 

Ajj = diag(l/x D ) AjjdiagCxp ) 

A 

• - dlagd/^ ) ^diagCXp ) 

1 1 
A 

B d -= diag(l/x D ) B^iagCup ) 

A 

B w = diagCl/x^) B^diagCup ) 

3. Calculate: 

“d - (1 - a)u D 1 /(u D j - U W j > 

6 J ' (aU D J - U W J ) ''S j - u W J ) 

4. Calculate: 


j " 1* . . . ,m 


1 "—i * 

k j ■ ( Vv 


id 


k l 3 - 


5. Calculate: 


m 


j 411 

•J- l 

j=i 


2 R 2m+1 


k 2m+l 1 " 


V 


C 2mf 2 V /x n 

i i 


1 r , 1 

S 2 - 2. kj_ 

2 j-l ^ 


i 1 1 
S 4 " k 2m+2 




j B 1, . . . )1D 

i ® 1 f • • • y n 


i c 1| • • # }ii 


6. Go to Algorithm II. 
Send: 


i i i i 
s l» s 2 , s 3 , s 4 , a, e 


Receive: r*. r*, r*. r J, y 


1 D lj • • ■ yii 


if-m 

t - 

i 


7. Calculate: 


'2m+l 


1 1 
C 2m+2 “ r 4 


3 K 1 ^ ■ jin 


■nRwirw- -fi— n-F-W - 
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.1 i 

-Ori-j i . C j 

i, 2 in 

V a -» 


kt - 


1 1 
r 3 C m+j 


1 * 1 f • • • ftl 


8* Output: 


c 2m+2 


V s j 

°i W i 

A A 

V \ 


j * 1, . . . ,tn 

1 ® 1 | • • • yti 


Algorithm II . 

1 . Input: s^,S2,S2>s^,e,a 

2. Calculate: 


Pi " s i “ 


S 3~ S 4 


1 "1 1-a 


Po “ s o “ 


V S 4 


2 2 1-a 


3. Minimize by line search: 


ax 


P 2 * P 1 


x-1 _ a -1 
a-1 


x - 


x i 

a -1 
a-1 


for -10 _< x <_ 10, x ^ 0, xi* 1 


4. Calculate: 


r 2 * 3? 


r 3 * 


8 „ - 


.'*-1 

r > ; 2 ~ 

s^, + r 3 (a -a) 

1— a ^ 


« ( V S 2 + r 2 r 3 (a 2 - 1 ” 
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- as 3 - r 3 (a -a) 
— 


5. Return to Algorithm 1.6 


3.4 


Straight Linear Affine Model 


As a check that the Power law term has significant effect. on the function 
g(u), a straight linear affine approximation to g(u) was generated. This 
model is then subjected to the same analysis as models 3A and 4. 


i 

3.5 Numerical Results 


[ 2 ] 


The algorithm of the previous section was applied to data obtained 


using DYNGEN with and u D specified as in Section 2. An off-design 

point was obtained using u^ - (.72727, .72727), with the resulting norm- 


alized state 

a 

• 

9000, .7897, .7381, 

.9401, . 

.9454). 

The normalized A 

and B matrices 

are 




— 




9“ 









-3.8 

-1.277 

2.067 

-1.152 

1.448 


-.00259 

.3553 

A 

2.748 

-5.39 

1.585 

-1.991 

1.071 

A 

.2116 

-.31618 

*D _ 

377.9 

49.51 

-264.9 

86.807 

78.91 

B w = 

12.54- 

-13.774 

31.26 

139.39 

-6.269 

-88.69 

27.83 

-.6201 

-99.3 


-176.5 

23.91 

-10.27 

-37.4 

-246.7 


157.78 

6.84 


— 







—4 




(3.5-1) 


-4.744 

-1.3888 

3.2468 -1.4591 

1.1969 


-.04546 

.0013 

.82.86 

-26.726 

2.5585 -1.8609 

.45548 

A 

.0086 

-.0121 

475.73 

137.55 

-328.91 27.791 

91.495 

V' 

2.434 

-.613 

-50.103 

110.91 

63.188 -116.69 

8.2883 

.67865 

-97.467 

-186.77 

-67.682 

-41.681 24.586 

-243.23 


203.44 

.64755 

— 





— 

— 


(3.5-2) 


Using the parameter value a = 
equilibrium functions &(u)* as 


.7, the Cj coefficients which specify the 
in Section 3 are given by the matrix 



.24267 

-.00218 

1.90082 

8.09916 

.02864 

.73088" 



1.01593 

.85407 

.89872 

.66919 - 

.81879 

-.05121 


c ■ 

. 73445 

.10133 

6.90586 

3.09409 

.011495 

.15272 

(3.5-3) 


.77234 

-.35905 

2.45867 

2.87415 - 

075198 

.66191 


_. 39503 

-.27262 

-3.44682 

13.4468 

.01838 

.85921 


This matrix together with the values 

A 

<*=1 . 1 and 

8=0.1 

and the matrices 


Ad and A w completely specify Model 3A. 


Another model which we will call Model 3B is easily obtained by using 
a li n e a r a ffine approximation to g(u) such that gfcy = * D> g«y » 

Moiel 3B is specified by a - e' 1 , n . 2.31778, 6 = -1.31778 and the co- 


efficient matrix 


.1553 .0028 1.0 1.0 0. .8418 

.1619 .1707 1.0 1.0 0. ..6674 

.5351 -.1208 1.0 1.0 0. .5857 

.5878 -.49313 1.0 1.0 0. .9053 

.2962 -.2099 1.0 1.0 0. .9137 


(3.5-4) 


CHAPTER IV 


QUASI-HERMITE INTERPOLATION MODEL 4 

4.1 Introduction 

A known method for matching a polynomial to the values and derivatives 
of a function at several points is Hermite interpolation. However, this 
method is formulated in general only for the one dimensional case. [3] Some 
works exist which apply this method to an n dimensional case, but only 
under certain narrow restrictions. [5] The single variable case, its restricted 
application in n dimensions, and a modified application in two dimensions, 
are discussed in this chapter. 

4.2 Hermite Interpolation for a Single Variable 

This presentation of the Hermite interpolation method is drawn from 
Hildebrand. [3] His notation is preserved, as closely-.as possible, here 
and in the resultant computer program. 

12 tn 

If the values of g(u) are known at m points, u = u , u , ...,u , 
define: 

tt(u) ** (u - u^)(u - u^) . . . (u - u m ) (4.2.1) 

and: 

^(u) (4.2.2) 

(U-U ) TT^U 1 ) 

with the properties: 

ir(u j ) =0 j = l,...,ra (4.2.3) 

and : 

^(u^) - i * l,...,m j = 1 m (4.2.4) 


183 
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where 6^ is the Kronecker delta (6^ = 1, 6^ - 0 for all i i 2 J). 

With these defined, the polynomial of degree m-1 which takes on the 

values g(u^) , g(u^) , . . . ,g(u m ) can be expressed as: 

m , , 

y(u) * l t (u)g(u K ) (4.2.5) 

k*=l 

Suppose both g(u) and g'(u) are known for u = u\ u^,...,u m , it 
is possible to determine a polynomial of degree 2m-l with these values 
and derivatives. We shall assume this polynomial is expressible in the 
form: 

y(u) = l h k (u) g(u k ) + l h k (u) g'(u k ) (4.2.6) 

k=l k=l 

i ~i 

where h (u) and h (u) , i = l,...,m are polynomials of maximum degree 
2m-l. The requirement y(u^) = g(u^) will be satisfied if: 

h*(u^) = 6 and h 1 ^) = 0 - (4.2.7) 

and the requirement y"(u^) = g'(u^) will be satisfied if: 

h '* 1 (u^ ) = 0 and h' 1 ^) = 6 (4.2.8) 

Since £ X (u) is a polynomial of degree m-1 which satisfied (4.2.4), 
i 2 

then [Z (u) ] is a polynomial of degree 2m-2 which satisfies (4.2.4) 
and whose derivative is zero at when i $ j . So if h*(u) and h i (u) 

are polynomials of degree 2m-l, then: 

h*(u) » r^uM-e^u)] 2 and h^u) = s^u) [^(u) ] 2 (4.2.9) 

i i 

where r (u) and s (u) are linear functions of u, so that (4.2.7) and 
(4.2.8) will be satisfied when i ^ j. These four conditions, when i = j, 
then yield: 
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r 1 (u^) * 1 r'V) +2 ^(u 1 ) - 0 
s*(u*) » 0 s^Cu 1 ) *» 1 

from which follows: 


(4.2.10) 

(4.2.11) 


r*(u) = 1-2 -d^(u*) (u “ u 1 ) and s*(u) «= (u - u*) (4.2.12) 

So; by combining (4.2.6), (4.2.9), and (4.2.12) we obtain the desired poly- 
nomial in the form: 


where: 


and: 


• ? k k I? -k k 

y(u) - l h (u) g(u) + l h (u) g"(u ) 
k=l k»l 

h*(u) = r*(u) f£*(u) and h*(u) = s^(u) [£*(u) 


(4.2.13) 


(4.2.14) 


r*(u) = 1-2 -£"*(u*)(u - u*) and s*(u) = (u - u^) (4.2.15) 

This result is known as Hermite' s interpolation formula, or the formula for 
osculating interpolation. 


4.3 Problems of Hermite Interpolation in n Dimensions 

If Hermite Interpolation were to be applied in n dimensions, the 
task would be to determine m sets of n+1 polynomials with properties sim- 
ilar to h and h. Specifically, if we assume that the desired polynomial 
can be expressed in the form: 

y(u) = l h k (u) g(u k ) + l l h^ k (u) (u k ) (4.3.1) 

k=l J-l k=l du J 

then these polynomials must have the properties: 

hV) = 6^ and h lk (u j ) = 0 . (4.3.2) 


and: 
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(Sh 1 , k. _ . 6h Jk , k. . 

• (u ) = 0 and (u ) = 


(4.3.3) 


<5u 


<5u. 


5 ’"J 

corresponding to the conditions (4.2.7) and (4.2.8). However, the further 
, condition: 


(u ) ■ 0 for all i ^ k 

du. 


(4.3.4) 


must also be satisfied. This final condition cannot be satisfied by the 
polynomials described in the previous section. In [5], JJiijima treats a 

i 

special case, in which the existence of certain orthogonal polynomials al- 
lows the application of Hermite interpolation to carefully chosen data in 
two dimensions. However, this method is not universally applicable to ar- 
bitrary data. 


4.4 Quasi-Hermite Approach for TWo Controls 

Given that no general method of Hermite interpolation in two variables 
was found, the following adaptation of the one dimensional case was applied. 
The value of control u 2 (nozzle area) was held constant at the design 
point value, and Hermite interpolation was applied to a set of points gen- 
erated by varying u^ (fuel flow). Both values, and derivatives with re- 
spect to u^ were matched at these data points. Then, for each value of 
u^, a value A was chosen, and control u 2 was varied by this amount, 
both plus and minus. A function was then chosen to match values at these 
new points, without altering the function at the original points. The re- 
sulting polynomial is of the form: 

y(u) * l h k (u 1 ) g(u k ,u*) 9 k (u 2 ) + l bNu^) g'(u k ,u 2 ) (4.4.1) 

k=l k“l 


v* 

... — - w -~ - 
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where: 


e (u 2 ) = 1 + [ 


, k 1 , .k. > k 1 .k N 

g(u^» u 2 + A ) - g^, u 2 - A ) 


2 A k g(u k ,u 2 > 


g(u k ,u 2 + A k ) + g(u k ,u 2 - A k ) - 2g(u k ,u 2 ) ^ 2 

+ f k .2 , k l, 

2(A ) gCu^u^) 


(4.4.2) 


This function has the property that: 


and: 


e k (u 2 ) = 1 when u 2 * u* 


0 k (u 2 ) = g(u k ,u 2 + A k )/g(u k ) when u 2 - U* + A k 
0 k (u 2 ) = g(u k ,u 2 - A k )/g(u k ) when u 2 = u* - A k 


(4.4.3) 


(4.4.4) 


and since h^(u^) = 5^, the resultant polynomial will match values and de- 
rivatives with respect to u^ along the u 2 = u 2 line, and values at u 2 = 

1 . A k ~~ 

u 2 ± A . 

4.5 Formation of A(x,u) 

Having chosen an approximation to g(u), it remains to choose a method 
for interpolating values for A(x,u) to complete the model. Lagrangian 
interpolation was used to match values only at three points along the u 2 - 
u 2 line, and at u = (uj\ u P + A"*") . The results of this approach are em- 
bodied in the following equations. First define: 


A1 

** A(x,u) 

at 

u = 

1 

u = 

, i 

<v 

lv 

u 2 ), 

x = 

gCu 1 ) 

A3 

= A(x,u) 

at 

u = 

3 

u = 

, 3 
(u l* 

3. 

V’ 

x - 

g(u 3 ) 

A5 

= A(x,u) 

at 

u = 

5 

u = 

( V 

5. 

U 2 } ’ 

X = 

g(u 5 ) 

AP 

* A(x,u) 

at 

u - 

P 

u = 

( V 

1 j. 
U 2 + 

A P ) ,x = 

g(u ? ) 

AM 

= A(x,.u) 

at 

u = 

M 

u = 


1 

U 2 “ 

A 1 ) ,x = 

g(u M ) 


(4.5.1) 
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( 4 . 5 . 2 ) 


( 4 . 5 . 3 ) 



CHAPTER V 


MODEL RESPONSE COMPARISONS 

5.1 Introduction 

Before subjecting the models to the Dynamic Programming Algorithm, 
some effort was made to examine their closeness of fit to DYNGEN data. 

Steady state values of models 1 and 4 are compared, and fuel flow step re- 
sponses of models 3A, 3B, and 4 are plotted in comparison to DYNGEN responses. 
A description of the step response program is also included. 

5.2 Steady State Comparison of Models 1 and 4 

The function g(u) represents a mapping from the control space into 
the state space, relating fixed controls to steady states. It is not only 
useful in the model form: 

x = A(x,u) [x - g(u)] (5.2-1) 

but should also approximate the operating line of the plant. 

Such a comparison is made here between g(u) for model 4 'and the 
DYNGEN simulator. Nozzle area was held constant, as fuel flow was varied 
from 9.0 to 1.1 by 0.02. All values are normalized . Percentage error is 
also computed, and shows the model's excellent agreement in its range of 
accuracy, and rapid deterioration outside that range. 
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Table 5.2-1 
X(l) - NC 


I 


j fuel flow 

DYNGEN 

Model 4 

% error 

0.90 

.97275 

.97288 

0.01 

r 0.92 

.97761 

.97790 

0.03 

i 0.94 

.98326 

.98326 

0.00 

0.96 

.98887 

.98892 

0.01 

0.98 

.99445 

.99461 

0.02 

1 1*00 

1.0000 

1.0000 

0.00 

** 1.02 

1.0051 

1.0051 

0.00 

1.04 

1.0102 

1.0113 

0.11 

1 1.06 

1.0152 

1.0230 

0.77 

1.08 

1.0201 

1.0513 

3.06 

1.10 

1.0244 

1.1195 

9.28 


I 

I 

I 

I 

I 

I 6 - -S 

I 

I 


Table 5.2-2 


X(2) * NF 


fuel flow 


0.90 

0.92 

0.94 

0.96 

1.00 

1.02 

1.04 

1.06 

1.08 

1.10 


DYNGEN 


.97132 

.97883 

.98427 

.98961 

1.0000 

1.0046 

1.0091 

1.0136 

1.0180 

1.0221 


Model 4 


.97099 

.97817 

.98425 

.98948 

1.0000 

1.0011 

.98046 

.89094 

.62787 

-.013230 


% error 


-0.03 

-0.07 

0.00 

- 0.01 

0.00 

-0.35 

-2.84 

- 12.10 

-38.32 

-101.29 
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Table 5.2-3 


X(3) - P4 


fuel flow 


0.90 

0.92 

0.94 

0.96 

0.98 

1.00 

1.02 

1.04 

1.06 

1.08 

1.10 


DYNGEN Model 4 


.92038 

.93623 

.95225 

.96821 

.98413 

1.0000 

1.0148 

1.0295 

1.0441 

1.0587 

1.0727 


.92042 

.93633 

.95225 

,96824 

;98434 

1.0000 

1.0125 

1.0136 

.98374 

.88135 

.62822 


% error 


0.00 

0.01 

0.00 

0.00 

0.02 

0.00 

-0.23 

-1.54 

-5,78 

-16.75 

-41.44 






Table 5.2-4 
X(4) - P7 
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Table 5.2-5 


X(5) *= U4 


lel flow 

DYNGEN 

0.90 

.96625 

0.92 

.97304 

0.94 

.97992 

0.96 

.98670 

0.98 

.99339 

1.00 

1.0000 

1.02 

1.0074 

1.04 

1.0147 

1.06 

1.0219 

1.08 

1.0290 

1.10 

1.0365 


Model 4 

X error 

.96624 

0.00 

.97308 

0.00 

.97992 

0.00 

.98665 

-0.01 

.99320 

-0.02 

1.0000 

0.00 

1.0089 

0.15 

1.0248 

1.00 

1.0573 

3.46 

1.1231 

9.14 

1.2494 

20.54 
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5 . 3 Program Layout 

The method chosen for generating time response data was a Euler in- 
egration with a user varied time step. After specifying initial controls, 
the user provides a control sequence of time step, duration (in iterations), 
and controls.. This structure allows the user to provide smaller time in- 
crements for the steeper portions of the response, and to alter the con- 
trols during the response. The step response program creates a file of 
time-state n-tuples., which are plotted against similar DYNGEN data by an- 
other program. 

5.4 Fuel Step Response for Models 3A, 3B, and 4 

Each of the three models was subjected to a fuel flow step from 0.8 
to 1.0, and the response plotted against the same response by DYNGEN. 

These graphs show that all three models match DYNGEN closely, but that 
model 4 is a better fit than either 3A or 3B. 
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Figure 5.4-2 
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